library(data.table)
library(truncnorm)
library(coda)
devtools::install_github("stan-dev/loo")
library(loo)
library(rstan) 
library(shinystan)
library(invgamma)
library(bayesplot)
library(posterior)

devtools::install_github("jfrench/bayesutils")

setwd('/Users/AM/Documents/_CU Masters/2020 fall Bayesian_7393/code/Bayesian_Statistics_Class_Code/Exam_code')
#df1 = load("diamonds_simple.rda")
#rm(df1, diamonds_simple)
df = bayesutils::diamonds_simple

df$lprice = log(df$price)
df$lcarat = log(df$carat)

lq_theta_y = function(sigmaSQ, beta0, beta1, lpr = df$lprice, lcar = df$lcarat) {
  ld = dnorm(beta0, 0, 100, log = TRUE) + 
    dnorm(beta1, 0, 100, log = TRUE) + 
    dinvgamma(sigmaSQ, shape = 0.01, rate = 0.01, log = TRUE) +
    sum(dnorm(lpr, mean = (beta0 + beta1 * lcar), sd = sqrt(sigmaSQ), log = TRUE))
  return(ld)
}

mh = function(B, theta_start) {
  theta = array(0, c((B+1), 3), dimnames = list(c(), c("sigmaSQ", "beta0", "beta1")))
  
  theta[1,1] = theta_start[1]
  theta[1,2] = theta_start[2]
  theta[1,3] = theta_start[3]

  for (i in 2:dim(theta)[1]) {
    
    ###  step for sigmaSQ
    beta0_star = theta[(i-1),2]
    beta1_star = theta[(i-1),3]
    
    sigmaSQ_star = rtruncnorm(n = 1, a=0, b=Inf, mean = theta[(i-1),1], sd = 0.1)

    num_logr = lq_theta_y(sigmaSQ = sigmaSQ_star, beta0 = beta0_star, beta1 = beta1_star) -
      log(dtruncnorm(x = sigmaSQ_star, a=0, b=Inf, mean = theta[(i-1),1], sd = 0.1))
    den_logr = lq_theta_y(sigmaSQ = theta[i-1, 1], beta0 = beta0_star, beta1 = beta1_star) -
      log(dtruncnorm(x = theta[i-1, 1], a=0, b=Inf, mean = sigmaSQ_star, sd = 0.1)) 
## !!! is the mean correct? Should  we use in den_logr mean = sigmaSQ_star or mean = theta[i-1, 1]? 
# check c-componentwise-mh v1 line 68, 72
    logr = num_logr - den_logr
    if (log(runif(1)) <= min(logr, 0)) {
      theta[i,1] = sigmaSQ_star
    } else {
      theta[i,1] = theta[(i - 1), 1]
    }
    
    ###  step for beta0
    beta1_star = theta[(i-1),3] # it is the repeated code, but I need it to keep the interpretability
    sigmaSQ_star = theta[i,1] # update sigmaSQ_star after the Gibbs step for sigmaSQ
    
    beta0_star = rnorm(1, theta[(i-1),2], 0.1) # !!! check the parametrization (0.1 or 0.1^2)
    
    num_logr = lq_theta_y(sigmaSQ = sigmaSQ_star, beta0 = beta0_star, beta1 = beta1_star) -
      dnorm(x = beta0_star, mean = theta[(i-1),2], sd = 0.1, log = TRUE)
    den_logr = lq_theta_y(sigmaSQ = sigmaSQ_star, beta0 = theta[i-1, 2], beta1 = beta1_star) -
      dnorm(x = theta[(i-1),2], mean = beta0_star, sd = 0.1, log = TRUE)
## !!! is the mean correct? Should  we use in den_logr mean = sigmaSQ_star or mean = theta[i-1, 1]? 
# check c-componentwise-mh v1 line 68, 72
    logr = num_logr - den_logr
    if (log(runif(1)) <= min(logr, 0)) {
      theta[i,2] = beta0_star
    } else {
      theta[i,2] = theta[(i - 1), 2]
    }
    
    ###  step for beta1
    sigmaSQ_star = theta[i,1] # it is the repeated code, but I need it to keep the interpretability
    beta0_star = theta[i,2] # update beta0_star after the Gibbs step for beta0

    beta1_star = rnorm(1, theta[(i-1),3], 0.1) # !!! check the parametrization (0.1 or 0.1^2)
    
    num_logr = lq_theta_y(sigmaSQ = sigmaSQ_star, beta0 = beta0_star, beta1 = beta1_star) -
      dnorm(x = beta1_star, mean = theta[(i-1),3], sd = 0.1, log = TRUE)
    den_logr = lq_theta_y(sigmaSQ = sigmaSQ_star, beta0 = beta0_star, beta1 = theta[i-1, 3]) -
      dnorm(x = theta[(i-1),3], mean = beta1_star, sd = 0.1, log = TRUE)
## !!! is the mean correct? Should  we use in den_logr mean = sigmaSQ_star or mean = theta[i-1, 1]? 
# check c-componentwise-mh v1 line 68, 72
    logr = num_logr - den_logr
    if (log(runif(1)) <= min(logr, 0)) {
      theta[i,3] = beta1_star
    } else {
      theta[i,3] = theta[(i - 1), 3]
    }
    
  }
  return(theta)
}

B = 10^5 
keep = (B/2 + 1):(B + 1)
chain1 = mh(B, theta_start = c(0.1, -1, -1))
chain2 = mh(B, theta_start = c(0.3, 0, 0))
chain3 = mh(B, theta_start = c(0.5, -1, 1))
chain4 = mh(B, theta_start = c(0.2, 1, 1))

mc = mcmc.list(mcmc(chain1[keep,]), mcmc(chain2[keep,]),
                 mcmc(chain3[keep,]), mcmc(chain4[keep,]))
summary(mc)

Iterations = 1:50001
Thinning interval = 1 
Number of chains = 4 
Sample size per chain = 50001 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

          Mean      SD  Naive SE Time-series SE
sigmaSQ 0.1152 0.00925 2.068e-05      6.977e-05
beta0   8.3798 0.02240 5.008e-05      1.647e-04
beta1   1.5064 0.03135 7.011e-05      2.141e-04

2. Quantiles for each variable:

           2.5%    25%    50%    75%  97.5%
sigmaSQ 0.09853 0.1088 0.1145 0.1211 0.1348
beta0   8.33546 8.3649 8.3797 8.3948 8.4234
beta1   1.44522 1.4854 1.5064 1.5276 1.5680
gelman.diag(mc, autoburnin = FALSE)
Potential scale reduction factors:

        Point est. Upper C.I.
sigmaSQ       1.03       1.08
beta0         1.00       1.01
beta1         1.01       1.03

Multivariate psrf

1.03
# extracting the MCMC samples from the soda_fit object
samples = extract(r_fit_A)
ncycles = length(samples[[1]])
T <- length(df$price)

# each row of yrep is a sample from the pp distribution
yrep = matrix(0, ncol = T, nrow = ncycles)
for (i in seq_len(T)) {
  mui = samples$beta0 + samples$beta1 * df$lcarat[i]
  yrep[, i] = rnorm(ncycles, mean = mui, sd = sqrt(samples$sigmaSQ))
}

# approximate posterior predictive density for an observation
# with the same covariates as observation 1
#plot(density(yrep[,1]))
color_scheme_set("red")

# posterior predictive check
y = df$lprice
ppc_hist(y, yrep[sample(1:ncycles, 8),]) + ggtitle("Model A: histogram comparing the response to 8 replicated data sets")



# scatterplot of y vs average yrep
ppc_scatter_avg(y, yrep) + ggtitle("Model A: scatterplot comparing the response to the average replicated data set")


ppc_dens_overlay(y, yrep = yrep[sample(1:ncycles, 100),]) + ggtitle("Model A: density plot comparing the density of the response to the densities of 100 replicated data sets")


ppc_error_scatter_avg(y, yrep = yrep) + ggtitle("Model A: scatter plot comparing the response to the average predictive error")


ppc_intervals(y, yrep = yrep)


ppc_error_scatter_avg_vs_x(y = y, yrep = yrep, x = df$lcarat)

ppc_intervals(y = y, yrep = yrep, x = df$lcarat) + ggtitle("Model A: Posterior predictive intervals to the observed data values") + xlab("log carat") + ylab("log price")

# extracting the MCMC samples from the soda_fit object
samples = extract(r_fit_B)
ncycles = length(samples[[1]])
T <- length(df$price)

# each row of yrep is a sample from the pp distribution
yrep = matrix(0, ncol = T, nrow = ncycles)
for (i in seq_len(T)) {
  mui = samples$beta0 + samples$beta1 * df$carat[i] + samples$beta2 * df$carat[i] * df$carat[i]
  yrep[, i] = rnorm(ncycles, mean = mui, sd = sqrt(samples$sigmaSQ))
}

# approximate posterior predictive density for an observation
# with the same covariates as observation 1
#plot(density(yrep[,1]))
color_scheme_set("green")

# posterior predictive check
y = df$price
ppc_hist(y, yrep[sample(1:ncycles, 8),]) + ggtitle("Model B: histogram comparing the response to 8 replicated data sets")



# scatterplot of y vs average yrep
ppc_scatter_avg(y, yrep) + ggtitle("Model B: scatterplot comparing the response to the average replicated data set")


ppc_dens_overlay(y, yrep = yrep[sample(1:ncycles, 100),]) + ggtitle("Model B: density plot comparing the density of the response to the densities of 100 replicated data sets")


ppc_error_scatter_avg(y, yrep = yrep) + ggtitle("Model B: scatter plot comparing the response to the average predictive error")


ppc_intervals(y = y, yrep = yrep, x = df$carat) + ggtitle("Model B: Posterior predictive intervals to the observed data values") + xlab("carat") + ylab("price")

LS0tCnRpdGxlOiAiNzM5MyBFeGFtIDIgUiBDb2RlIEFuZHJlaSBNYXR2ZWV2IgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCmBgYHtyIERhdGEgYW5kIHNldCB1cH0KCmxpYnJhcnkoZGF0YS50YWJsZSkKbGlicmFyeSh0cnVuY25vcm0pCmxpYnJhcnkoY29kYSkKZGV2dG9vbHM6Omluc3RhbGxfZ2l0aHViKCJzdGFuLWRldi9sb28iKQpsaWJyYXJ5KGxvbykKbGlicmFyeShyc3RhbikgCmxpYnJhcnkoc2hpbnlzdGFuKQpsaWJyYXJ5KGludmdhbW1hKQpsaWJyYXJ5KGJheWVzcGxvdCkKbGlicmFyeShwb3N0ZXJpb3IpCgpkZXZ0b29sczo6aW5zdGFsbF9naXRodWIoImpmcmVuY2gvYmF5ZXN1dGlscyIpCgpzZXR3ZCgnL1VzZXJzL0FNL0RvY3VtZW50cy9fQ1UgTWFzdGVycy8yMDIwIGZhbGwgQmF5ZXNpYW5fNzM5My9jb2RlL0JheWVzaWFuX1N0YXRpc3RpY3NfQ2xhc3NfQ29kZS9FeGFtX2NvZGUnKQojZGYxID0gbG9hZCgiZGlhbW9uZHNfc2ltcGxlLnJkYSIpCiNybShkZjEsIGRpYW1vbmRzX3NpbXBsZSkKZGYgPSBiYXllc3V0aWxzOjpkaWFtb25kc19zaW1wbGUKCmRmJGxwcmljZSA9IGxvZyhkZiRwcmljZSkKZGYkbGNhcmF0ID0gbG9nKGRmJGNhcmF0KQoKYGBgCgpgYGB7ciBQcm9ibGVtIDEgbG9nIGRlbnNpdHkgRnVuY3Rpb259CgpscV90aGV0YV95ID0gZnVuY3Rpb24oc2lnbWFTUSwgYmV0YTAsIGJldGExLCBscHIgPSBkZiRscHJpY2UsIGxjYXIgPSBkZiRsY2FyYXQpIHsKICBsZCA9IGRub3JtKGJldGEwLCAwLCAxMDAsIGxvZyA9IFRSVUUpICsgCiAgICBkbm9ybShiZXRhMSwgMCwgMTAwLCBsb2cgPSBUUlVFKSArIAogICAgZGludmdhbW1hKHNpZ21hU1EsIHNoYXBlID0gMC4wMSwgcmF0ZSA9IDAuMDEsIGxvZyA9IFRSVUUpICsKICAgIHN1bShkbm9ybShscHIsIG1lYW4gPSAoYmV0YTAgKyBiZXRhMSAqIGxjYXIpLCBzZCA9IHNxcnQoc2lnbWFTUSksIGxvZyA9IFRSVUUpKQogIHJldHVybihsZCkKfQoKYGBgCgpgYGB7ciBQcm9ibGVtIDEgbWggRnVuY3Rpb259CgptaCA9IGZ1bmN0aW9uKEIsIHRoZXRhX3N0YXJ0KSB7CiAgdGhldGEgPSBhcnJheSgwLCBjKChCKzEpLCAzKSwgZGltbmFtZXMgPSBsaXN0KGMoKSwgYygic2lnbWFTUSIsICJiZXRhMCIsICJiZXRhMSIpKSkKICAKICB0aGV0YVsxLDFdID0gdGhldGFfc3RhcnRbMV0KICB0aGV0YVsxLDJdID0gdGhldGFfc3RhcnRbMl0KICB0aGV0YVsxLDNdID0gdGhldGFfc3RhcnRbM10KCiAgZm9yIChpIGluIDI6ZGltKHRoZXRhKVsxXSkgewogICAgCiAgICAjIyMgIHN0ZXAgZm9yIHNpZ21hU1EKICAgIGJldGEwX3N0YXIgPSB0aGV0YVsoaS0xKSwyXQogICAgYmV0YTFfc3RhciA9IHRoZXRhWyhpLTEpLDNdCiAgICAKICAgIHNpZ21hU1Ffc3RhciA9IHJ0cnVuY25vcm0obiA9IDEsIGE9MCwgYj1JbmYsIG1lYW4gPSB0aGV0YVsoaS0xKSwxXSwgc2QgPSAwLjEpCgogICAgbnVtX2xvZ3IgPSBscV90aGV0YV95KHNpZ21hU1EgPSBzaWdtYVNRX3N0YXIsIGJldGEwID0gYmV0YTBfc3RhciwgYmV0YTEgPSBiZXRhMV9zdGFyKSAtCiAgICAgIGxvZyhkdHJ1bmNub3JtKHggPSBzaWdtYVNRX3N0YXIsIGE9MCwgYj1JbmYsIG1lYW4gPSB0aGV0YVsoaS0xKSwxXSwgc2QgPSAwLjEpKQogICAgZGVuX2xvZ3IgPSBscV90aGV0YV95KHNpZ21hU1EgPSB0aGV0YVtpLTEsIDFdLCBiZXRhMCA9IGJldGEwX3N0YXIsIGJldGExID0gYmV0YTFfc3RhcikgLQogICAgICBsb2coZHRydW5jbm9ybSh4ID0gdGhldGFbaS0xLCAxXSwgYT0wLCBiPUluZiwgbWVhbiA9IHNpZ21hU1Ffc3Rhciwgc2QgPSAwLjEpKSAKIyMgISEhIGlzIHRoZSBtZWFuIGNvcnJlY3Q/IFNob3VsZCAgd2UgdXNlIGluIGRlbl9sb2dyIG1lYW4gPSBzaWdtYVNRX3N0YXIgb3IgbWVhbiA9IHRoZXRhW2ktMSwgMV0/IAojIGNoZWNrIGMtY29tcG9uZW50d2lzZS1taCB2MSBsaW5lIDY4LCA3MgogICAgbG9nciA9IG51bV9sb2dyIC0gZGVuX2xvZ3IKICAgIGlmIChsb2cocnVuaWYoMSkpIDw9IG1pbihsb2dyLCAwKSkgewogICAgICB0aGV0YVtpLDFdID0gc2lnbWFTUV9zdGFyCiAgICB9IGVsc2UgewogICAgICB0aGV0YVtpLDFdID0gdGhldGFbKGkgLSAxKSwgMV0KICAgIH0KICAgIAogICAgIyMjICBzdGVwIGZvciBiZXRhMAogICAgYmV0YTFfc3RhciA9IHRoZXRhWyhpLTEpLDNdICMgaXQgaXMgdGhlIHJlcGVhdGVkIGNvZGUsIGJ1dCBJIG5lZWQgaXQgdG8ga2VlcCB0aGUgaW50ZXJwcmV0YWJpbGl0eQogICAgc2lnbWFTUV9zdGFyID0gdGhldGFbaSwxXSAjIHVwZGF0ZSBzaWdtYVNRX3N0YXIgYWZ0ZXIgdGhlIEdpYmJzIHN0ZXAgZm9yIHNpZ21hU1EKICAgIAogICAgYmV0YTBfc3RhciA9IHJub3JtKDEsIHRoZXRhWyhpLTEpLDJdLCAwLjEpICMgISEhIGNoZWNrIHRoZSBwYXJhbWV0cml6YXRpb24gKDAuMSBvciAwLjFeMikKICAgIAogICAgbnVtX2xvZ3IgPSBscV90aGV0YV95KHNpZ21hU1EgPSBzaWdtYVNRX3N0YXIsIGJldGEwID0gYmV0YTBfc3RhciwgYmV0YTEgPSBiZXRhMV9zdGFyKSAtCiAgICAgIGRub3JtKHggPSBiZXRhMF9zdGFyLCBtZWFuID0gdGhldGFbKGktMSksMl0sIHNkID0gMC4xLCBsb2cgPSBUUlVFKQogICAgZGVuX2xvZ3IgPSBscV90aGV0YV95KHNpZ21hU1EgPSBzaWdtYVNRX3N0YXIsIGJldGEwID0gdGhldGFbaS0xLCAyXSwgYmV0YTEgPSBiZXRhMV9zdGFyKSAtCiAgICAgIGRub3JtKHggPSB0aGV0YVsoaS0xKSwyXSwgbWVhbiA9IGJldGEwX3N0YXIsIHNkID0gMC4xLCBsb2cgPSBUUlVFKQojIyAhISEgaXMgdGhlIG1lYW4gY29ycmVjdD8gU2hvdWxkICB3ZSB1c2UgaW4gZGVuX2xvZ3IgbWVhbiA9IHNpZ21hU1Ffc3RhciBvciBtZWFuID0gdGhldGFbaS0xLCAxXT8gCiMgY2hlY2sgYy1jb21wb25lbnR3aXNlLW1oIHYxIGxpbmUgNjgsIDcyCiAgICBsb2dyID0gbnVtX2xvZ3IgLSBkZW5fbG9ncgogICAgaWYgKGxvZyhydW5pZigxKSkgPD0gbWluKGxvZ3IsIDApKSB7CiAgICAgIHRoZXRhW2ksMl0gPSBiZXRhMF9zdGFyCiAgICB9IGVsc2UgewogICAgICB0aGV0YVtpLDJdID0gdGhldGFbKGkgLSAxKSwgMl0KICAgIH0KICAgIAogICAgIyMjICBzdGVwIGZvciBiZXRhMQogICAgc2lnbWFTUV9zdGFyID0gdGhldGFbaSwxXSAjIGl0IGlzIHRoZSByZXBlYXRlZCBjb2RlLCBidXQgSSBuZWVkIGl0IHRvIGtlZXAgdGhlIGludGVycHJldGFiaWxpdHkKICAgIGJldGEwX3N0YXIgPSB0aGV0YVtpLDJdICMgdXBkYXRlIGJldGEwX3N0YXIgYWZ0ZXIgdGhlIEdpYmJzIHN0ZXAgZm9yIGJldGEwCgogICAgYmV0YTFfc3RhciA9IHJub3JtKDEsIHRoZXRhWyhpLTEpLDNdLCAwLjEpICMgISEhIGNoZWNrIHRoZSBwYXJhbWV0cml6YXRpb24gKDAuMSBvciAwLjFeMikKICAgIAogICAgbnVtX2xvZ3IgPSBscV90aGV0YV95KHNpZ21hU1EgPSBzaWdtYVNRX3N0YXIsIGJldGEwID0gYmV0YTBfc3RhciwgYmV0YTEgPSBiZXRhMV9zdGFyKSAtCiAgICAgIGRub3JtKHggPSBiZXRhMV9zdGFyLCBtZWFuID0gdGhldGFbKGktMSksM10sIHNkID0gMC4xLCBsb2cgPSBUUlVFKQogICAgZGVuX2xvZ3IgPSBscV90aGV0YV95KHNpZ21hU1EgPSBzaWdtYVNRX3N0YXIsIGJldGEwID0gYmV0YTBfc3RhciwgYmV0YTEgPSB0aGV0YVtpLTEsIDNdKSAtCiAgICAgIGRub3JtKHggPSB0aGV0YVsoaS0xKSwzXSwgbWVhbiA9IGJldGExX3N0YXIsIHNkID0gMC4xLCBsb2cgPSBUUlVFKQojIyAhISEgaXMgdGhlIG1lYW4gY29ycmVjdD8gU2hvdWxkICB3ZSB1c2UgaW4gZGVuX2xvZ3IgbWVhbiA9IHNpZ21hU1Ffc3RhciBvciBtZWFuID0gdGhldGFbaS0xLCAxXT8gCiMgY2hlY2sgYy1jb21wb25lbnR3aXNlLW1oIHYxIGxpbmUgNjgsIDcyCiAgICBsb2dyID0gbnVtX2xvZ3IgLSBkZW5fbG9ncgogICAgaWYgKGxvZyhydW5pZigxKSkgPD0gbWluKGxvZ3IsIDApKSB7CiAgICAgIHRoZXRhW2ksM10gPSBiZXRhMV9zdGFyCiAgICB9IGVsc2UgewogICAgICB0aGV0YVtpLDNdID0gdGhldGFbKGkgLSAxKSwgM10KICAgIH0KICAgIAogIH0KICByZXR1cm4odGhldGEpCn0KCmBgYAoKYGBge3IgUHJvYmxlbSAxIG1oIFJ1bn0KCkIgPSAxMF41IAprZWVwID0gKEIvMiArIDEpOihCICsgMSkKY2hhaW4xID0gbWgoQiwgdGhldGFfc3RhcnQgPSBjKDAuMSwgLTEsIC0xKSkKY2hhaW4yID0gbWgoQiwgdGhldGFfc3RhcnQgPSBjKDAuMywgMCwgMCkpCmNoYWluMyA9IG1oKEIsIHRoZXRhX3N0YXJ0ID0gYygwLjUsIC0xLCAxKSkKY2hhaW40ID0gbWgoQiwgdGhldGFfc3RhcnQgPSBjKDAuMiwgMSwgMSkpCgptYyA9IG1jbWMubGlzdChtY21jKGNoYWluMVtrZWVwLF0pLCBtY21jKGNoYWluMltrZWVwLF0pLAogICAgICAgICAgICAgICAgIG1jbWMoY2hhaW4zW2tlZXAsXSksIG1jbWMoY2hhaW40W2tlZXAsXSkpCnN1bW1hcnkobWMpCgpgYGAKCmBgYHtyIFByb2JsZW0gMiB0byBhc3Nlc3MgdGhlIGNvbnZlcmdlbmNlfQoKa2VlcCA9IChCLzIgKyAxOTAwMSk6KEIvMiArIDIwMDAxKQoKbWMgPSBtY21jLmxpc3QobWNtYyhjaGFpbjFba2VlcCxdKSwgbWNtYyhjaGFpbjJba2VlcCxdKSwKICAgICAgICAgICAgICAgICBtY21jKGNoYWluM1trZWVwLF0pLCBtY21jKGNoYWluNFtrZWVwLF0pKQpjb2RhOjp0cmFjZXBsb3QobWMpCmNvZGE6OmF1dG9jb3JyLnBsb3QobWMsIGxhZy5tYXggPSAxMDAsIGF1dG8ubGF5b3V0ID0gVFJVRSkKZ2VsbWFuLmRpYWcobWMsIGF1dG9idXJuaW4gPSBGQUxTRSkKZ2V3ZWtlLmRpYWcobWMpCgpgYGAKCmBgYHtyIFByb2JsZW0gMyBNb2RlbCBBLCBpbmNsdWRlPUZBTFNFfQpzdGFuX21vZF9BID0gIgpkYXRhIHsKICBpbnQgVDsKICB2ZWN0b3JbVF0gbHByOyAgICAvLyBsb2cgcHJpY2UgCiAgdmVjdG9yW1RdIGxjYXI7ICAgLy8gbG9nIGNhcmF0IAogIHJlYWw8bG93ZXI9MD4gdjsgIC8vIHNhbXBsZSB2YXJpYW5jZSBvZiBsb2cgcHJpY2UgCgp9CgpwYXJhbWV0ZXJzIHsKICByZWFsPGxvd2VyPTA+IHNpZ21hU1E7CiAgcmVhbCBiZXRhMDsKICByZWFsIGJldGExOwp9Cgp0cmFuc2Zvcm1lZCBwYXJhbWV0ZXJzIHsKICB2ZWN0b3JbVF0gbXU7CiAgZm9yIChpIGluIDE6VCkgewogICAgbXVbaV0gPSBiZXRhMCArIGJldGExICogbGNhcltpXTsKICB9Cgp9Cgptb2RlbCB7CiAgc2lnbWFTUSB+IGludl9nYW1tYSgwLjAxLCAwLjAxKTsKICBiZXRhMCB+IG5vcm1hbCgwLCAxMDApOwogIGJldGExIH4gbm9ybWFsKDAsIDEwMCk7CgogIGZvciAoaSBpbiAxOlQpCiAgICBscHJbaV0gfiBub3JtYWwobXVbaV0sIHNxcnQoc2lnbWFTUSkpOwp9CmdlbmVyYXRlZCBxdWFudGl0aWVzIHsKICB2ZWN0b3JbVF0gbG9nX2xpazsKICB2ZWN0b3JbVF0geV9yZXA7CgogIC8vdmVjdG9yW1RdIGxpazsKICByZWFsIFJic3E7ICAgICAgICAgICAgICAvLyBnb29kbmVzcy1vZi1maXQKICBSYnNxID0gMSAtIHNpZ21hU1EvdjsKICAKICBmb3IgKGkgaW4gMTpUKSB7CiAgICBsb2dfbGlrW2ldID0gbm9ybWFsX2xwZGYobHByW2ldIHwgbXVbaV0sIHNxcnQoc2lnbWFTUSkpOwogICAgeV9yZXBbaV0gPSBub3JtYWxfcm5nKG11W2ldLCBzcXJ0KHNpZ21hU1EpKTsKICAgIC8vbGlrW2ldID0gZXhwKGxvZ19saWtbaV0pOwogIH0KCn0KIgoKVCA8LSBsZW5ndGgoZGYkbHByaWNlKQp2ID0gdmFyKGRmJGxwcmljZSkKc2V0LnNlZWQoOTApCgptb2RlbF9uYW1lID0gIk1vZGVsX0FfIgpzdGFuX2RhdF9BID0gbGlzdChUID0gVCwgbHByID0gZGYkbHByaWNlLCBsY2FyID0gZGYkbGNhcmF0LCB2ID0gdikKcl9maXRfQSA9IHN0YW4obW9kZWxfY29kZSA9IHN0YW5fbW9kX0EsIGRhdGEgPSBzdGFuX2RhdF9BLCAKICAgICAgICAgICAgIGl0ZXIgPSA1MDAwLCBjaGFpbnMgPSA0LCAKICAgICAgICAgICAgIGNvbnRyb2wgPSBsaXN0KG1heF90cmVlZGVwdGggPSAyMSkpCgojZmlsZV9uYW1lID0gcGFzdGUobW9kZWxfbmFtZSwgYXMuY2hhcmFjdGVyKFN5cy5EYXRlKCkpLCAiLnJkYSIsIHNlcD0iIikKI3NldHdkKCcvVXNlcnMvQU0vRG9jdW1lbnRzL19DVSBNYXN0ZXJzLzIwMjAgZmFsbCBCYXllc2lhbl83MzkzL2NvZGUvQmF5ZXNpYW5fU3RhdGlzdGljc19DbGFzc19Db2RlL0V4YW1fY29kZScpCiNyZWFkUkRTKHJfZml0LCBmaWxlID0gIiIpIAojc2F2ZVJEUyhyX2ZpdCwgZmlsZSA9IGZpbGVfbmFtZSwgY29tcHJlc3MgPSAieHoiKSAKCnN1bW1hcnkocl9maXRfQSwgcGFycyA9IGMoInNpZ21hU1EiLCAiYmV0YTAiLCAiYmV0YTEiKSwgcHJvYiA9IGMoMC4wMjUsIDAuOTc1KSkkc3VtbWFyeQoKI3NzbyA8LSBsYXVuY2hfc2hpbnlzdGFuKHJfZml0X0EpCgpgYGAKCmBgYHtyIFByb2JsZW0gMyBNb2RlbCBCLCBpbmNsdWRlPUZBTFNFfQpzdGFuX21vZF9CID0gIgpkYXRhIHsKICBpbnQgVDsKICB2ZWN0b3JbVF0gcHI7ICAgICAvLyBwcmljZSAKICB2ZWN0b3JbVF0gY2FyOyAgICAvLyBjYXJhdCAKICByZWFsPGxvd2VyPTA+IHY7ICAvLyBzYW1wbGUgdmFyaWFuY2Ugb2YgbG9nIHByaWNlIAoKfQoKcGFyYW1ldGVycyB7CiAgcmVhbDxsb3dlcj0wPiBzaWdtYVNROwogIHJlYWwgYmV0YTA7CiAgcmVhbCBiZXRhMTsKICByZWFsIGJldGEyOwp9Cgp0cmFuc2Zvcm1lZCBwYXJhbWV0ZXJzIHsKICB2ZWN0b3JbVF0gbXU7CiAgZm9yIChpIGluIDE6VCkgewogICAgbXVbaV0gPSBiZXRhMCArIGJldGExICogY2FyW2ldICsgYmV0YTIgKiBjYXJbaV0gKiBjYXJbaV07CiAgfQoKfQoKbW9kZWwgewogIHNpZ21hU1EgfiBpbnZfZ2FtbWEoMC4wMSwgMC4wMSk7CiAgYmV0YTAgfiBub3JtYWwoMCwgMTAwKTsKICBiZXRhMSB+IG5vcm1hbCgwLCAxMDApOwogIGJldGEyIH4gbm9ybWFsKDAsIDEwMCk7CgogIGZvciAoaSBpbiAxOlQpCiAgICBwcltpXSB+IG5vcm1hbChtdVtpXSwgc3FydChzaWdtYVNRKSk7Cn0KCmdlbmVyYXRlZCBxdWFudGl0aWVzIHsKICB2ZWN0b3JbVF0geV9yZXA7CiAgdmVjdG9yW1RdIGxvZ19saWs7CiAgLy92ZWN0b3JbVF0gbGlrOwogIHJlYWwgUmJzcTsgICAgICAgICAgICAgIC8vIGdvb2RuZXNzLW9mLWZpdAogIFJic3EgPSAxIC0gc2lnbWFTUS92OwogIAogIGZvciAoaSBpbiAxOlQpIHsKICAgIGxvZ19saWtbaV0gPSBub3JtYWxfbHBkZihwcltpXSB8IG11W2ldLCBzcXJ0KHNpZ21hU1EpKTsKICAgIHlfcmVwW2ldID0gbm9ybWFsX3JuZyhtdVtpXSwgc3FydChzaWdtYVNRKSk7CgogICAgLy9saWtbaV0gPSBleHAobG9nX2xpa1tpXSk7CiAgfQoKfQoiCgpUIDwtIGxlbmd0aChkZiRwcmljZSkKdiA9IHZhcihkZiRwcmljZSkKc2V0LnNlZWQoOTEpCgptb2RlbF9uYW1lID0gIk1vZGVsX0JfIgpzdGFuX2RhdF9CID0gbGlzdChUID0gVCwgcHIgPSBkZiRwcmljZSwgY2FyID0gZGYkY2FyYXQsIHYgPSB2KQpyX2ZpdF9CID0gc3Rhbihtb2RlbF9jb2RlID0gc3Rhbl9tb2RfQiwgZGF0YSA9IHN0YW5fZGF0X0IsIAogICAgICAgICAgICAgaXRlciA9IDUwMDAsIGNoYWlucyA9IDQsIAogICAgICAgICAgICAgY29udHJvbCA9IGxpc3QobWF4X3RyZWVkZXB0aCA9IDIxKSkKCiNmaWxlX25hbWUgPSBwYXN0ZShtb2RlbF9uYW1lLCBhcy5jaGFyYWN0ZXIoU3lzLnRpbWUoKSksICIucmRhIiwgc2VwPSIiKQojc2V0d2QoJy9Vc2Vycy9BTS9Eb2N1bWVudHMvX0NVIE1hc3RlcnMvMjAyMCBmYWxsIEJheWVzaWFuXzczOTMvY29kZS9CYXllc2lhbl9TdGF0aXN0aWNzX0NsYXNzX0NvZGUvRXhhbV9jb2RlJykKI3NhdmVSRFMocl9maXQsIGZpbGUgPSBmaWxlX25hbWUsIGNvbXByZXNzID0gInh6IikgCnN1bW1hcnkocl9maXRfQiwgcGFycyA9IGMoInNpZ21hU1EiLCAiYmV0YTAiLCAiYmV0YTEiLCAiYmV0YTIiKSwgcHJvYiA9IGMoMC4wMjUsIDAuOTc1KSkkc3VtbWFyeQoKc3NvIDwtIGxhdW5jaF9zaGlueXN0YW4ocl9maXRfQikKCmBgYAoKYGBge3IgUHJvYmxlbSA0IHBwX2NoZWNrIE1vZGVsIEF9CiMgZXh0cmFjdGluZyB0aGUgTUNNQyBzYW1wbGVzIGZyb20gdGhlIHNvZGFfZml0IG9iamVjdApzYW1wbGVzID0gZXh0cmFjdChyX2ZpdF9BKQpuY3ljbGVzID0gbGVuZ3RoKHNhbXBsZXNbWzFdXSkKVCA8LSBsZW5ndGgoZGYkcHJpY2UpCgojIGVhY2ggcm93IG9mIHlyZXAgaXMgYSBzYW1wbGUgZnJvbSB0aGUgcHAgZGlzdHJpYnV0aW9uCnlyZXAgPSBtYXRyaXgoMCwgbmNvbCA9IFQsIG5yb3cgPSBuY3ljbGVzKQpmb3IgKGkgaW4gc2VxX2xlbihUKSkgewogIG11aSA9IHNhbXBsZXMkYmV0YTAgKyBzYW1wbGVzJGJldGExICogZGYkbGNhcmF0W2ldCiAgeXJlcFssIGldID0gcm5vcm0obmN5Y2xlcywgbWVhbiA9IG11aSwgc2QgPSBzcXJ0KHNhbXBsZXMkc2lnbWFTUSkpCn0KCiMgYXBwcm94aW1hdGUgcG9zdGVyaW9yIHByZWRpY3RpdmUgZGVuc2l0eSBmb3IgYW4gb2JzZXJ2YXRpb24KIyB3aXRoIHRoZSBzYW1lIGNvdmFyaWF0ZXMgYXMgb2JzZXJ2YXRpb24gMQojcGxvdChkZW5zaXR5KHlyZXBbLDFdKSkKY29sb3Jfc2NoZW1lX3NldCgicmVkIikKCiMgcG9zdGVyaW9yIHByZWRpY3RpdmUgY2hlY2sKeSA9IGRmJGxwcmljZQpwcGNfaGlzdCh5LCB5cmVwW3NhbXBsZSgxOm5jeWNsZXMsIDgpLF0pICsgZ2d0aXRsZSgiTW9kZWwgQTogaGlzdG9ncmFtIGNvbXBhcmluZyB0aGUgcmVzcG9uc2UgdG8gOCByZXBsaWNhdGVkIGRhdGEgc2V0cyIpCgoKIyBzY2F0dGVycGxvdCBvZiB5IHZzIGF2ZXJhZ2UgeXJlcApwcGNfc2NhdHRlcl9hdmcoeSwgeXJlcCkgKyBnZ3RpdGxlKCJNb2RlbCBBOiBzY2F0dGVycGxvdCBjb21wYXJpbmcgdGhlIHJlc3BvbnNlIHRvIHRoZSBhdmVyYWdlIHJlcGxpY2F0ZWQgZGF0YSBzZXQiKQoKcHBjX2RlbnNfb3ZlcmxheSh5LCB5cmVwID0geXJlcFtzYW1wbGUoMTpuY3ljbGVzLCAxMDApLF0pICsgZ2d0aXRsZSgiTW9kZWwgQTogZGVuc2l0eSBwbG90IGNvbXBhcmluZyB0aGUgZGVuc2l0eSBvZiB0aGUgcmVzcG9uc2UgdG8gdGhlIGRlbnNpdGllcyBvZiAxMDAgcmVwbGljYXRlZCBkYXRhIHNldHMiKQoKcHBjX2Vycm9yX3NjYXR0ZXJfYXZnKHksIHlyZXAgPSB5cmVwKSArIGdndGl0bGUoIk1vZGVsIEE6IHNjYXR0ZXIgcGxvdCBjb21wYXJpbmcgdGhlIHJlc3BvbnNlIHRvIHRoZSBhdmVyYWdlIHByZWRpY3RpdmUgZXJyb3IiKQoKcHBjX2ludGVydmFscyh5ID0geSwgeXJlcCA9IHlyZXAsIHggPSBkZiRsY2FyYXQpICsgZ2d0aXRsZSgiTW9kZWwgQTogUG9zdGVyaW9yIHByZWRpY3RpdmUgaW50ZXJ2YWxzIHRvIHRoZSBvYnNlcnZlZCBkYXRhIHZhbHVlcyIpICsgeGxhYigibG9nIGNhcmF0IikgKyB5bGFiKCJsb2cgcHJpY2UiKQpgYGAKCmBgYHtyIFByb2JsZW0gNCBwcF9jaGVjayBNb2RlbCBCfQojIGV4dHJhY3RpbmcgdGhlIE1DTUMgc2FtcGxlcyBmcm9tIHRoZSBzb2RhX2ZpdCBvYmplY3QKc2FtcGxlcyA9IGV4dHJhY3Qocl9maXRfQikKbmN5Y2xlcyA9IGxlbmd0aChzYW1wbGVzW1sxXV0pClQgPC0gbGVuZ3RoKGRmJHByaWNlKQoKIyBlYWNoIHJvdyBvZiB5cmVwIGlzIGEgc2FtcGxlIGZyb20gdGhlIHBwIGRpc3RyaWJ1dGlvbgp5cmVwID0gbWF0cml4KDAsIG5jb2wgPSBULCBucm93ID0gbmN5Y2xlcykKZm9yIChpIGluIHNlcV9sZW4oVCkpIHsKICBtdWkgPSBzYW1wbGVzJGJldGEwICsgc2FtcGxlcyRiZXRhMSAqIGRmJGNhcmF0W2ldICsgc2FtcGxlcyRiZXRhMiAqIGRmJGNhcmF0W2ldICogZGYkY2FyYXRbaV0KICB5cmVwWywgaV0gPSBybm9ybShuY3ljbGVzLCBtZWFuID0gbXVpLCBzZCA9IHNxcnQoc2FtcGxlcyRzaWdtYVNRKSkKfQoKIyBhcHByb3hpbWF0ZSBwb3N0ZXJpb3IgcHJlZGljdGl2ZSBkZW5zaXR5IGZvciBhbiBvYnNlcnZhdGlvbgojIHdpdGggdGhlIHNhbWUgY292YXJpYXRlcyBhcyBvYnNlcnZhdGlvbiAxCiNwbG90KGRlbnNpdHkoeXJlcFssMV0pKQpjb2xvcl9zY2hlbWVfc2V0KCJncmVlbiIpCgojIHBvc3RlcmlvciBwcmVkaWN0aXZlIGNoZWNrCnkgPSBkZiRwcmljZQpwcGNfaGlzdCh5LCB5cmVwW3NhbXBsZSgxOm5jeWNsZXMsIDgpLF0pICsgZ2d0aXRsZSgiTW9kZWwgQjogaGlzdG9ncmFtIGNvbXBhcmluZyB0aGUgcmVzcG9uc2UgdG8gOCByZXBsaWNhdGVkIGRhdGEgc2V0cyIpCgoKIyBzY2F0dGVycGxvdCBvZiB5IHZzIGF2ZXJhZ2UgeXJlcApwcGNfc2NhdHRlcl9hdmcoeSwgeXJlcCkgKyBnZ3RpdGxlKCJNb2RlbCBCOiBzY2F0dGVycGxvdCBjb21wYXJpbmcgdGhlIHJlc3BvbnNlIHRvIHRoZSBhdmVyYWdlIHJlcGxpY2F0ZWQgZGF0YSBzZXQiKQoKcHBjX2RlbnNfb3ZlcmxheSh5LCB5cmVwID0geXJlcFtzYW1wbGUoMTpuY3ljbGVzLCAxMDApLF0pICsgZ2d0aXRsZSgiTW9kZWwgQjogZGVuc2l0eSBwbG90IGNvbXBhcmluZyB0aGUgZGVuc2l0eSBvZiB0aGUgcmVzcG9uc2UgdG8gdGhlIGRlbnNpdGllcyBvZiAxMDAgcmVwbGljYXRlZCBkYXRhIHNldHMiKQoKcHBjX2Vycm9yX3NjYXR0ZXJfYXZnKHksIHlyZXAgPSB5cmVwKSArIGdndGl0bGUoIk1vZGVsIEI6IHNjYXR0ZXIgcGxvdCBjb21wYXJpbmcgdGhlIHJlc3BvbnNlIHRvIHRoZSBhdmVyYWdlIHByZWRpY3RpdmUgZXJyb3IiKQoKcHBjX2ludGVydmFscyh5ID0geSwgeXJlcCA9IHlyZXAsIHggPSBkZiRjYXJhdCkgKyBnZ3RpdGxlKCJNb2RlbCBCOiBQb3N0ZXJpb3IgcHJlZGljdGl2ZSBpbnRlcnZhbHMgdG8gdGhlIG9ic2VydmVkIGRhdGEgdmFsdWVzIikgKyB4bGFiKCJjYXJhdCIpICsgeWxhYigicHJpY2UiKQpgYGAKCmBgYHtyIFByb2JsZW0gNSBNb2RlbCBDLCBpbmNsdWRlPUZBTFNFfQpzdGFuX21vZF9DID0gIgpkYXRhIHsKICBpbnQgVDsKICB2ZWN0b3JbVF0gcHI7ICAgIC8vICBwcmljZSAKICB2ZWN0b3JbVF0gbGNhcjsgICAvLyBsb2cgY2FyYXQgCiAgcmVhbDxsb3dlcj0wPiB2OyAgLy8gc2FtcGxlIHZhcmlhbmNlIG9mIGxvZyBwcmljZSAKCn0KCnBhcmFtZXRlcnMgewogIHJlYWw8bG93ZXI9MD4gc2lnbWFTUTsKICByZWFsIGJldGEwOwogIHJlYWwgYmV0YTE7Cn0KCnRyYW5zZm9ybWVkIHBhcmFtZXRlcnMgewogIHZlY3RvcltUXSBtdTsKICBmb3IgKGkgaW4gMTpUKSB7CiAgICBtdVtpXSA9IGJldGEwICsgYmV0YTEgKiBsY2FyW2ldOwogIH0KCn0KCm1vZGVsIHsKICBzaWdtYVNRIH4gaW52X2dhbW1hKDAuMDEsIDAuMDEpOwogIGJldGEwIH4gbm9ybWFsKDAsIDEwMCk7CiAgYmV0YTEgfiBub3JtYWwoMCwgMTAwKTsKCiAgZm9yIChpIGluIDE6VCkKICAgIHByW2ldIH4gbG9nbm9ybWFsKG11W2ldLCBzcXJ0KHNpZ21hU1EpKTsKfQpnZW5lcmF0ZWQgcXVhbnRpdGllcyB7CiAgdmVjdG9yW1RdIGxvZ19saWs7CiAgdmVjdG9yW1RdIHlfcmVwOwoKICAvL3ZlY3RvcltUXSBsaWs7CiAgcmVhbCBSYnNxOyAgICAgICAgICAgICAgLy8gZ29vZG5lc3Mtb2YtZml0CiAgUmJzcSA9IDEgLSBzaWdtYVNRL3Y7CiAgCiAgZm9yIChpIGluIDE6VCkgewogICAgbG9nX2xpa1tpXSA9IGxvZ25vcm1hbF9scGRmKHByW2ldIHwgbXVbaV0sIHNxcnQoc2lnbWFTUSkpOwogICAgeV9yZXBbaV0gPSBsb2dub3JtYWxfcm5nKG11W2ldLCBzcXJ0KHNpZ21hU1EpKTsKICAgIC8vbGlrW2ldID0gZXhwKGxvZ19saWtbaV0pOwogIH0KCn0KIgoKVCA8LSBsZW5ndGgoZGYkcHJpY2UpCnYgPSB2YXIoZGYkcHJpY2UpCnNldC5zZWVkKDkyKQoKbW9kZWxfbmFtZSA9ICJNb2RlbF9DXyIKc3Rhbl9kYXRfQyA9IGxpc3QoVCA9IFQsIHByID0gZGYkcHJpY2UsIGxjYXIgPSBkZiRsY2FyYXQsIHYgPSB2KQpyX2ZpdF9DID0gc3Rhbihtb2RlbF9jb2RlID0gc3Rhbl9tb2RfQywgZGF0YSA9IHN0YW5fZGF0X0MsIAogICAgICAgICAgICAgaXRlciA9IDUwMDAsIGNoYWlucyA9IDQsIAogICAgICAgICAgICAgY29udHJvbCA9IGxpc3QobWF4X3RyZWVkZXB0aCA9IDIxKSkKCiNmaWxlX25hbWUgPSBwYXN0ZShtb2RlbF9uYW1lLCBhcy5jaGFyYWN0ZXIoU3lzLkRhdGUoKSksICIucmRhIiwgc2VwPSIiKQojc2V0d2QoJy9Vc2Vycy9BTS9Eb2N1bWVudHMvX0NVIE1hc3RlcnMvMjAyMCBmYWxsIEJheWVzaWFuXzczOTMvY29kZS9CYXllc2lhbl9TdGF0aXN0aWNzX0NsYXNzX0NvZGUvRXhhbV9jb2RlJykKI3JlYWRSRFMocl9maXQsIGZpbGUgPSAiIikgCiNzYXZlUkRTKHJfZml0LCBmaWxlID0gZmlsZV9uYW1lLCBjb21wcmVzcyA9ICJ4eiIpIAoKc3VtbWFyeShyX2ZpdF9DLCBwYXJzID0gYygic2lnbWFTUSIsICJiZXRhMCIsICJiZXRhMSIpLCBwcm9iID0gYygwLjAyNSwgMC45NzUpKSRzdW1tYXJ5Cgojc3NvIDwtIGxhdW5jaF9zaGlueXN0YW4ocl9maXRfQSkKCmBgYAoKYGBge3IgUHJvYmxlbSA1IE1vZGVsIEQsIGluY2x1ZGU9RkFMU0V9CnN0YW5fbW9kX0QgPSAiCmRhdGEgewogIGludCBUOwogIHZlY3RvcltUXSBwcjsgICAgIC8vICBwcmljZSAKICB2ZWN0b3JbVF0gbGNhcjsgICAvLyBsb2cgY2FyYXQgCiAgdmVjdG9yW1RdIG5vdGljOyAgLy8gIG5vdGljZWFibGUgCgogIHJlYWw8bG93ZXI9MD4gdjsgIC8vIHNhbXBsZSB2YXJpYW5jZSBvZiBsb2cgcHJpY2UgCgp9CgpwYXJhbWV0ZXJzIHsKICByZWFsPGxvd2VyPTA+IHNpZ21hU1E7CiAgcmVhbCBiZXRhMDsKICByZWFsIGJldGExOwogIHJlYWwgYmV0YTI7Cn0KCnRyYW5zZm9ybWVkIHBhcmFtZXRlcnMgewogIHZlY3RvcltUXSBtdTsKICBmb3IgKGkgaW4gMTpUKSB7CiAgICBtdVtpXSA9IGJldGEwICsgYmV0YTEgKiBsY2FyW2ldICsgYmV0YTIgKiBub3RpY1tpXTsKICB9Cgp9Cgptb2RlbCB7CiAgc2lnbWFTUSB+IGludl9nYW1tYSgwLjAxLCAwLjAxKTsKICBiZXRhMCB+IG5vcm1hbCgwLCAxMDApOwogIGJldGExIH4gbm9ybWFsKDAsIDEwMCk7CiAgYmV0YTIgfiBub3JtYWwoMCwgMTAwKTsKCgogIGZvciAoaSBpbiAxOlQpCiAgICBwcltpXSB+IGxvZ25vcm1hbChtdVtpXSwgc3FydChzaWdtYVNRKSk7Cn0KZ2VuZXJhdGVkIHF1YW50aXRpZXMgewogIHZlY3RvcltUXSBsb2dfbGlrOwogIHZlY3RvcltUXSB5X3JlcDsKCiAgLy92ZWN0b3JbVF0gbGlrOwogIHJlYWwgUmJzcTsgICAgICAgICAgICAgIC8vIGdvb2RuZXNzLW9mLWZpdAogIFJic3EgPSAxIC0gc2lnbWFTUS92OwogIAogIGZvciAoaSBpbiAxOlQpIHsKICAgIGxvZ19saWtbaV0gPSBsb2dub3JtYWxfbHBkZihwcltpXSB8IG11W2ldLCBzcXJ0KHNpZ21hU1EpKTsKICAgIHlfcmVwW2ldID0gbG9nbm9ybWFsX3JuZyhtdVtpXSwgc3FydChzaWdtYVNRKSk7CiAgICAvL2xpa1tpXSA9IGV4cChsb2dfbGlrW2ldKTsKICB9Cgp9CiIKClQgPC0gbGVuZ3RoKGRmJHByaWNlKQp2ID0gdmFyKGRmJHByaWNlKQpzZXQuc2VlZCg5MykKCm1vZGVsX25hbWUgPSAiTW9kZWxfRF8iCnN0YW5fZGF0X0QgPSBsaXN0KFQgPSBULCBwciA9IGRmJHByaWNlLCBsY2FyID0gZGYkbGNhcmF0LCAKICAgICAgICAgICAgICAgICAgbm90aWMgPSBhcy5pbnRlZ2VyKGRmJGluY2x1c2lvbnMgPT0gIm5vdGljZWFibGUiKSAsIHYgPSB2KQpyX2ZpdF9EID0gc3Rhbihtb2RlbF9jb2RlID0gc3Rhbl9tb2RfRCwgZGF0YSA9IHN0YW5fZGF0X0QsIAogICAgICAgICAgICAgaXRlciA9IDUwMDAsIGNoYWlucyA9IDQsIAogICAgICAgICAgICAgY29udHJvbCA9IGxpc3QobWF4X3RyZWVkZXB0aCA9IDIxKSkKCiNmaWxlX25hbWUgPSBwYXN0ZShtb2RlbF9uYW1lLCBhcy5jaGFyYWN0ZXIoU3lzLkRhdGUoKSksICIucmRhIiwgc2VwPSIiKQojc2V0d2QoJy9Vc2Vycy9BTS9Eb2N1bWVudHMvX0NVIE1hc3RlcnMvMjAyMCBmYWxsIEJheWVzaWFuXzczOTMvY29kZS9CYXllc2lhbl9TdGF0aXN0aWNzX0NsYXNzX0NvZGUvRXhhbV9jb2RlJykKI3JlYWRSRFMocl9maXQsIGZpbGUgPSAiIikgCiNzYXZlUkRTKHJfZml0LCBmaWxlID0gZmlsZV9uYW1lLCBjb21wcmVzcyA9ICJ4eiIpIAoKc3VtbWFyeShyX2ZpdF9ELCBwYXJzID0gYygic2lnbWFTUSIsICJiZXRhMCIsICJiZXRhMSIsICJiZXRhMiIpLCBwcm9iID0gYygwLjAyNSwgMC45NzUpKSRzdW1tYXJ5Cgojc3NvIDwtIGxhdW5jaF9zaGlueXN0YW4ocl9maXRfQSkKCmBgYAoKYGBge3IgUHJvYmxlbSA2LCBpbmNsdWRlPUZBTFNFfQp3YWljX2xvbyA9IGZ1bmN0aW9uIChmaXQpIHsKICBsb2dfbGlrID0gZXh0cmFjdF9sb2dfbGlrKGZpdCwgbWVyZ2VfY2hhaW5zID0gRkFMU0UpCiAgcl9lZmYgPSBleHAocmVsYXRpdmVfZWZmKGxvZ19saWspKQogIHdsID0gYyh3YWljKGxvZ19saWspJHdhaWMsIGxvbyhsb2dfbGlrLCByX2VmZiA9IHJfZWZmKSRsb29pYykKICByZXR1cm4od2wpCn0Kd2wgPSBhcnJheSgwLCBkaW09YygyLDMpLCAKICAgICAgICAgICBkaW1uYW1lcyA9IGxpc3QoYygiV0FJQyIsICJMT09JQyIpLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgYygiTW9kZWwgQiIsIk1vZGVsIEMiLCJNb2RlbCBEIikpKQp3bFsxLDFdID0gd2FpY19sb28ocl9maXRfQilbMV0Kd2xbMiwxXSA9IHdhaWNfbG9vKHJfZml0X0IpWzJdCndsWzEsMl0gPSB3YWljX2xvbyhyX2ZpdF9DKVsxXQp3bFsyLDJdID0gd2FpY19sb28ocl9maXRfQylbMl0Kd2xbMSwzXSA9IHdhaWNfbG9vKHJfZml0X0QpWzFdCndsWzIsM10gPSB3YWljX2xvbyhyX2ZpdF9EKVsyXQoKd2wKYGBgCgpgYGB7cn0KCmBgYAoK